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ABSTRACT: 

In this report, service routines are developed for linear, quadratic 
and cubic finite elements for two and three dimensional regions. Also, an 
equation solver of very large capacity is developed for systems of linear 
equations where the matrix of coefficients is symmetrical, positive def- 
inite and tightly banded along the min diagonal. 
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INTRODUCTION 



The stress-strain analysis of solid propellant rocket motors has 
been the object of numerous investigations in recent years. Various 
approximate methods of analysis have been proposed. The Finite Ele- 
ment method, without a doubt, showed the greatest promise of providing 
the information needed in the evaluation of real motors . In principle 
the method can provide analyses of any three dimensional bodies. Any 
number of material and geometrical discontinuities can be admitted. 

The method also allows for almost any type of static, dynamic and 
thermal loadings to be applied. In practice, however, this has proved 
to be a formidable computational task even for the largest digital com- 
puters in existence. 

An elegant compromise was obtained with the development of axi- 

★ 

symmetric ring elements (1) . However this compromise eliminated the 
possibility of solving for real motors with star shaped grains. In those 
cases it has generally been assumed that a plane strain analysis using 
elements of triangular shape would provide enough information for de- 
sign. The constant strain finite element triangle has been used ad 
nauseam in these studies, in spite of the fact that many other elements 
have been shown to provide much better results for the same computational 
effort (2). When the perforation changes shape, the combination of ring 

* 

Numbers in parenthesis refer to the reference list at the end. 
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elements and constant strain triangle cannot give reliable results. 

Recently the introduction of isoparametric elements in the literature 
has given a new impetus to attempts at solving real three-dimensional 
motors (3), (4), (5). It is to such a task that we address ourselves. 

In this work we shall not attempt to account for real viscoelastic be- 
havior, what we shall present is a linear elastic analysis method for 
three-dimensional objects. We will show that isoparametric elements 
can be economically generated in the computer, even for the monstrous 
(96 x 96), 32 nodal point bricks. A massive system of linear equations 
with a large half-bandwidth normally results from the analysis of real- 
istic problems. For these systems of equations we propose a new 
solver (6) . 



2 



Chapter I. Isoparametric Elements 



In this chapter we describe isoparametric elements and the numeri- 
cal techniques used to construct them. 

1 . 1 Advantages of isoparametric elements 

The development of isoparametric elements marked an important 
breakthrough in the finite element method (4). Until then, a great 
variety of elements were developed to solve specific types of problems. 
The most successful elements satisfied all the requirements to insure 
convergence to an exact result by mesh refinement. However the solu- 
tion of many problems required a combination of various elements which 
could not satisfy all the required conditions for successful convergence, 
when assembled to solve one problem. An example of such a situation 
is: an axisymmetric system where truncated conical ring elements are 
used to model a case around a rocket motor and ring elements of tri- 
angular cross section are used to model the grain (1). The usual coni- 
cal element admits transverse displacements given by a cubic poly- 
nomial where as the grain ring elements will allow only linear dis- 
placement in the same direction. Obviously, adjacent ring and case 
elements cannot maintain continuity of displacements under strain. 

Isoparametric element families on the other hand can be developed 
to yield compatible elements of any type: line, surface or volume. 
Furthermore these elements admit curved boundaries, resulting in accu- 
rate modeling of complex shapes even with very few elements. The 
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results obtained in references (3), (4) and (5) show that good results 
are obtained with very coarse mesh. 

1 . 2 Stiffness matrix of a finite element 

The strain energy of an element having a stiffness matrix [K] can 
be written as: 

u s = \ <V CK] {u.} (1) 

where the vector {u.} represents all the nodal coordinates of the element. 

The same energy can also be written from fundamental principles 
as follows: 

U s = 2.[ <ff i >tf l) dv (2) 

V 

where {cr} is a vector of all the stress components and a vector 
of corresponding strain components in the volume V. The constitutive 
laws can be written in matrix form as follows: 



fCT.3 - [*]{€,} 



(3) 



where [6] is a symmetric matrix of material properties. The strain dis- 



placement relations are also written in matrix form as follows: 

/ \ 
u(x, y, z) 



{€.} = [<p] < v(x, y, z) 
w(x, y, z) 






(4) 



where matrix [<p] is in general a rectangular matrix of differential opera- 



tors, and the vector (u(x, y, z), v(x, y, z), w(x, y, z)> represents 
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the three components of displacements in an appropriate global system 
of reference, (x, y, z). Finally the displacement functions are as- 
sumed to be of the following form: 

/ \ 



u(x, y, z) 

< v(x, y, z) \ = 
w(x, y, z) 



[n] {u,} 



(5) 



\ / 

where [h] is a rectangular matrix of carefully selected functions of 
(x, y, z). These functions must be selected to insure compatibility 
across common boundaries, to insure that rigid body motions will result 
in strainless states and finally to admit at least constant straining 
states . 

After substitution of (5), (4) and (3) into (2) and using the definition: 



[n*] = [<p] [h] 



( 6 ) 



we obtain the following result: 

U S = 2 <V (J Cn*] T [eKn*] dv)[u.} (7) 

The stiffness matrix is then: 

[K] = f [h*] T [e][h*] dv (8) 

J v 

1 .3 The isoparametric concept 

In isoparametric elements, the displacements are obtained in terms 
of shape functions defined in a system of curvilinear coordinates particu- 
lar to one element. 
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( 9 ) 



u(x, y , z) = N (g jt i,e )u 
v(x, y, z) = N.(£,n,c)v. 
w(x, y, z) = N ( ? „ 5 )w 

1 5 5 1 

where N^( 5,n,s) are the shape functions and u^, v , are nodal 
values of the displacement components in the global system (x, y, z). 
The global system of reference is then related to the curvilinear co- 
ordinates by the same shape functions as follows: 

x ( E ,n ,?) = N iU,n,C ) x i 

y(5,n,c) = N.U,n,t; )y. (10) 

z (5,n,c) =N i(c,n,; ) z i 

where x. , y. and z. are the nodal coordinates of the element in the 
ii l 

global system (x, y, z). It is not our intention here to trace the history 
of such a concept. The readers are referred to references (3) and (4) 
for the details and justification of the concept. The remaining task is 
to construct the shape functions. We have taken these functions from 
reference (4). They will be found in section 1.6, 1.7 and 1.8. 

1.4 The Jacobian Matrix 

In equation (6) , the derivatives of the shape functions in terms of 
the coordinates of the global system (x, y, z) are needed. These opera- 
tions require the Jacobian of the transformation: 
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where the square matrix is the Jacobian [j]. Using equation (10) the 
Jacobian matrix is easily evaluated as follows: 



[J] = 



dN 3N 2 s 

TT ' "aT 








x i y i z i 


5N 3N 2 




X 2 y 2 Z 2 


3rj b V 
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SN 3N 2 
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3C ' ac ' ***;, 




... 



( 12 ) 



When the Jacobian is known the required derivatives are easily ob- 
tained. The element of volume is then transformed as follows: 



dV = det (J) d£ drj dC (13) 

At this time everything is ideally suited for Gaussian numerical inte- 
gration. This will be explained later. 

1 . 5 Shape functions 

As mentioned previously we take the displacement functions directly 
from reference (4) . We write these shape functions for three dimensional 
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elements since they automatically contain the information for two and 



one dimensional elements. We also restrict our attention to tri 
tri-quadratic and tri-cubic polynomial expressions 
1 . 6 Linear elements 

The shape functions are: 

N. = £ (l + 4 )(i + rj Ml + C ) i = l, 2 8 

1 0 0 0 0 

where 4 = ££. and £. is ± 1, similarly for r) and £ . 

0 1 1 0 0 



2 (- 1 , 1 , 1 ) 



« (-1,1 ,- I) 

8 (1,-1, -I) 5(1,1 ,-l) 

Figure 1. Linear element (8 nodal points) 




1 . 7 Quadratic elements 

See figure 2 for the identification of the nodal points. 
Corner nodes: 1, 3, 5, 7, 13, 15, 17 and 19 

n. (i + 4 )(i + t) )(i + c )(4 + t) +r - 

lo o o oooo 

where £ = ££. and £. is ± 1, similarly for r] and £ . 

o 1 1 oo 

Mid side nodes: 2, 6, 18 and 14 

N. = 7 (i - 4 2 )(1 + n Mi + £ ) 

1 H oo 



linear. 



(14) 



(15) 



(16) 
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Mid side nodes: 4, 8, 20 and 16 



N. =J (1 ' ^ 2)(1 + C)d + 4 ) 
14 oo 

Mid side nodes: 9, 10, 11 and 12 

N = j(l - C 2 )(l + 4 )(1 + rj ) 

1 4 00 



5 4 3 




Figure 2 . Quadratic element (20 nodal points) 

1 . 8 Cubic elements 

See figure 3 for the identification of the nodal points . 
Corner nodes 1, 4 , 7, 10, 21, 24, 27 and 30 

N. = 77 (1 + £ )(1 + rj )(1 + C >C9(4 2 + T) 2 + C 2 

1 oq 0 0 0 

Mid side nodes: 2, 9, 29, 22, 3, 8, 28 and 23 

N. = TT U - 4 2 )U + 94 )(1 + r? )(1 + c ) 

1 000 

where £ = ± ~ , r j. = ± 1 and C . = ± 1 

s i 3 '1 1 



( 17 ) 



(18) 



) - 19] (19) 



( 20 ) 
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Mid side nodes: 12 , 32 , 25, 5, 11 , 31, 26 and 6 



N. = -J7- (1 - T? 2 ) ( 1 + 9rj )(1 + C )(1 + % ) (21) 

1 ^ 0 0 0 

where 77 . = ±“",£.=±1 and C . = ± 1 

'1 3 1 1 

Mid side nodes: 13, 14, 15, 16, 17, 18, 19 and 20 

N. - (1 = C 2 )(l + 9? )(1 + t )(1 + r? ) (22) 

1 0 0 0 

where C. = ± ~, E. = ±1 and 77 . = ± 1 . 

1 3 1 '1 

7 6 5 4 




Figure 3. Cubic element (32 nodal points) 



In all of these formulas if one coordinate , for example C = C = 1 » 

0 

the formulas degenerate to the corresponding correct formulas for the 
two dimensional element in plane (£, 77 ). The same is true for line 
elements . 
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1 . 9 Numerical integration 



When all the information contained in sections 1.4 to 1.8 is 
substituted into equation (8) we finally obtain an expression of the 
following form: 



In this last expression the matrix k(£, n> £) * s composed of elements 
in which the inverse of the Jacobian matrix is involved, furthermore, 
all the elements must be multiplied by the determinant of the Jacobian 
matrix. Obviously the process of reduction to the form shown in 
equation (23) is not trivial, however, a computational algorithm is not 
too difficult to obtain. 

Once the form of equation (23) is obtained, Gaussian integration 
is straightforward. The number of Gauss points to be used in this 
integration is not easy to determine, especially if the element has 
irregular boundaries . For trirectangular elements , two, four and five 
Gauss ordinates in each direction of the local system of reference 
will give exact results for the linear, quadratic and cubic elements 
respectively. For elements having curved boundaries the same num- 
ber of Gauss ordinates can only give approximate values to the stiff- 
ness matrix. We experimented with various numbers of Gauss points 
in those cases where the results were approximate, we observed that 
the dominant eigenvalues of the stiffness matrices were not sensitive 



1 1 1 




(23) 
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to the number of Gauss ordinates. We also observed that the trace 



of an approximate stiffness matrix was always smaller than the trace 
of the exact matrix. By increasing the number of Gauss points, the 
trace of the approximate matrix always approached the exact value 
from the low side. This error may actually be beneficial, since the 
trace of a stiffness matrix is known to be always greater than the 
trace of the actual stiffness of a real structure. For this reason we 
have adopted to use 2, 4 and 5 Gauss ordinates in the generation of 
the stiffness matrices. It will be necessary to show with actual 
problems whether or not this compromise is valid in general. 
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Chapter II. Computer Programs for the Generation of Stiffness Matrices 



In this chapter we describe several computer programs used to 
evaluate the stiffness matrices of linear, quadratic and cubic ele- 
ments for two and three dimensional bodies. 

2 . 1 Stiffness matrices for two dimensional bodies 

Appendix one contains the listing of a computer program that was 
used to evaluate the stiffness matrices of two-dimensional elements. 
Although the listing contains only one quadrature formula with five 
Gauss ordinates, an exhaustive series of tests was conducted with 
formulas using from two to ten Gauss ordinates. The five point formula 
was finally adopted for practical considerations; the formation times 
for one element were: 0.3, 1.0, 2.5 seconds for linear, quadratic 
and cubic elements respectively, and this was judged to be adequate 
for this phase of our research. 

Following the strategy developed in section 1.9, two sub- 
routines are needed to do all the work. The first subroutine named 
QUAD 5 performs the numerical integration, the second named FORMK 
assembles the stiffness matrix k(£, rj, £) where (£, r\, C) take their 
values at one of the Gauss points. Upon exit from the subroutines, the 
user will obtain the stiffness matrix of one of the elements in this 
family and also the transformation matrices between nodal strains 
(stresses) and nodal displacements. The calling statement is given 
in the appendix. 
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2.2 Stiffness matrices for three dimensional bodies 



Appendix two contains three separate programs for the generation 
of the stiffnesses of linear, quadratic and cubic elements for three 
dimensional bodies. The three programs are similar in structure to 
the above program, they were coded separately for convenience during 
the testing phase of this study. The FORTRAN calling statements are 
all contained in the appendix. The number of Gauss point used was 
determined in section 1.9 for these elements. 

Comments distributed throughout these programs should make the 
codes intelligible to an experienced user of finite element codes. 

2 . 3 Numerical testing of stiffness matrices 

The subroutines mentioned in the previous two articles were 
tested under the control of another program not listed here. Eigen- 
values and eigenvectors were calculated for all these elements. Two 
dimensional elements always had three zero eigenvalues and three 
dimensional elements always had six zero eigenvalues. Static equili- 
brium checks were also applied to all elements and equilibrium was 
always satisfied. 

2 .4 Hybrid elements 

A remarkable property of isoparametric elements is that once a 
high order element is obtained it can easily be degenerated to a lower 
order by simple matrix transformations. For example let us consider 
the eight nodal point two-dimensional element: 
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(24) 



{F} = [K] {u.} 

where [K] is the stiffness matrix and {u^} is the vector of nodal dis- 
placements; 

(u.} = <u 1 .v 1 ,u 2 ,v 2 .u 3 ,v 3 ,u 4 ,u 4 ,u 5 ,v 5 ,u 6 ,v 6 ,u 7 ,v 7 ,u 8 ,v 8 >’ 

(25) 



see figure 4. 

To degenerate this element to the corresponding four nodal point 
element, it is only necessary to impose a simple constraint to the mid 
points of the element. If the mid-points are such that their nodal dis- 
placements are equal to the average of the displacements of their 
corresponding corner points the following transformation matrix is 
readily constructed: 

{u. } = [T] {u. } (26) 

★ iji 

where {u.] = (u^v^u^v^u^v^u^v^ (27) 



2 I 3 





— —HI 


► • — « 






<> ' 
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< > ■■■+ 
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7 5 



Figure 4, Hybrid elements 
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The matrix [T ] in this case is simply: 
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The stiffness matrix for the four nodal point element is simply: 

[K*] = [T] T [K] [T] (29) 

If a five nodal point element is desired, for example, suppose we 
would like to keep nodal point 8 in the above example, then two columns 
would be added to the matrix [T 3 - Column 9 of [ T] would have a one 
in row 15 and column 10 a one in row 16, all the other terms of rows 
15 and 16 would be zeros. The congruent transformation (29) would 
give the stiffness of the hybrid element. 

With such elements it is possible to combine different types of 
elements in the solution of one problem and still preserve all the 
conditions guaranteeing convergence to an exact result by mesh size 
reduction. 
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Chapter III. An Equation Solver of Very Large Capacity 



In this chapter a very simple equation solver of very large capa- 
city is developed. The method used does not require a very large in- 
core memory; however, a fast random access device like a magnetic 
disk drive is required. 

3 . 1 Introduction 

The solution of very large systems of linear equations, 

{R} = [K] * {r} (30) 

where [K] is symmetric and positive definite, is an activity that must 
be entered by the stress analyst when using the finite element method. 
When the first automatic stress analysis programs were developed, the 
Gauss-Seidel iteration method, with various over relaxation schemes,, 
was almost universally employed to solve the equations. It was soon 
realized that, because of the properties of [K], a direct Gaussian 
elimination algorithm would be much more economical of computer time. 
However, matrix [K], must be in core during elimination, and this puts 
a very stringent limitation on the total number of equations that can 
be solved. 

For the purpose of comparing various methods, we will presume 
that we have at our disposition a machine with 32,000 words of core 
storage. Of the total memory, let us assume that approximately 
25,000 words are available for the storage of [K]. For such a machine, 
then, approximately 160 equations, at most, can be solved by this 
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method . 



Because the matrix [K] is usually sparse, a special ordering of 
the unknowns can produce a matrix [K] that is tightly banded about 
the main diagonal. Taking advantage of this character of [K], several 
programs were devised that required only (N.M) words for the storage 
of [K]. (See Figure 5 for the meaning of N and M). The short program 
listed in Appendix 3 is of that character. If M is very small, then 
the number of equations that can be solved can become fairly large. 

For example, using the same machine as above, if M = 50, then 
N = 500. 

With the development of two and three dimensional elements, 
however, this improvement proved insufficient. These problems fre- 
quently lead to systems of several thousands of equations with a 
half band width in the hundreds . Equation solvers were then developed 
that did not require to have all the equations in core. Professor E. 

Wilson of the University of California at Berkeley (U . S.) several 
years ago developed a band solver that could be used with tape stor- 
age to solve large systems. The method used required (2M.M) words 
of in-core storage for [K] and the total number of equations that could 
be solved was only limited by the capacity of a tape. For the machine 
used above, the maximum value of M is approximately 110. Subsequently, 
as mentioned in Reference (8) , several workers developed algorithms 
requiring only (M + 2)(M + 3)/2 words of core storage. For the machine 
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used above, the maximum value of M is approximately 220. 

With large systems of equations, another difficulty becomes 
evident. Round off starts to play a very important role. An obvious 
remedy is to convert everything to double precision. However, this 
roughly doubles the core storage requirement. Round off can some- 
times be eliminated by using the same programs to calculate corrections 
to the answer by an iterative technique. However, such corrections 
are not always valid when the matrix [K] is locally ill-conditioned. 

For this reason, many workers prefer to use double precision all the 
way. The block solver presented here uses double precision every- 
where . 

As mentioned in References (9) and (10), no method has been 
found that is faster and more accurate than Gaussian elimination for 
a system in which the matrix of coefficients is dense. For dense 
system banded along the main diagonal, it is reasonable to conjecture 
that the same conclusion will hold if the algorithm is modified to 
eliminate all operations involving the off diagonal zeros. Improve- 
ments can only be accomplished by taking advantage of the sparsity 
within the band, as shown in References (11) and (12). These improve- 
ments, however, are obtained at the cost of a more elaborate code to 
recognize the particular type of sparsity that can lead to important 
reductions of computations. The economies of storage, on the other 
hand, are not sufficient to allow for the solution of very large systems. 
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Frontal solvers, as described in References (13) and (8), also produce 
savings of storage and execution time, but again prohibit the solution 
of very large systems. 

Our primary concern here was to design a solver that would permit 

the solution of very large systems even with a small computer. The 

block solver proposed in this paper does not pose any limits on the 

half band width nor on the total number of equations, as long as 

storage is available on a fast random access external device. The 

★ 

program contained in Appendix 4 requires 5800 bytes of storage for 
the object codes and a working space depending on the block size 
selected; for block sizes of (25, 50, 100) the storage requirements 
are (15600, 61200, 242400) bytes respectively. As mentioned above, 
a fast random access facility is required. In our computer center, 
we used IMB/2314 disk packs. On such a device, 55000 records, 

400 bytes long, can be stored. This is enough for 2,750,000 double 
words. For example, 5000 equations with a half band width of 451 
would fit with a complete load vector on one disk. 

M is the half band width 

N is the total number of 
equations 

Figure 5. Banded Character of K 
A byte is a word segment of 8 bits of information. 
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3 . 2 Elimination by blocks 



Consider the system of equations shown below: 
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(31) 



The first equation is solved for r^ to get: 



r i K 11 ^ R 1 K 12 r 2 K 13 r 3 “* _K lM r M^ 



After substitution in equations 2 to M, we obtain the following results: 



R 2 R 2 " K 12 K 11 R 1 



(33) 



R , = R.. - R. 

M M 1M 11 1 



(34) 



for the load vectors. All the coefficients in rows 2 to M from Column 2 
to M are reduced by similar operations; a typical example for term is: 
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( 35 ) 



K 34 - K 34 ' <3 K U K !4 
Then, Equation 2 is solved for and substituted in Equation 3 to M + 1 
and the entire process is repeated until we get only the last equation: 



R = K r 

N NN N 



(36) 



This last equation can now be solved for r^ and a simple process of 
back substitution in reverse order gives all the unknowns. 

The elimination process described above preserves the symmetry 
of [K]. This process is valid, whether the elements of [K] are 
individual coefficients or square submatrices of coefficients. From now 
on, we presume that the [K^] are submatrices that we will call blocks. 
A careful examination of the elimination steps reveals that no more 
than three different blocks of coefficients are needed. With this in 
mind, the rest of the task is trivial. The simple solver contained in 
Appendix 3 was modified to do the elimination by blocks. Appendix 4 
contains the listings of the four subroutines needed. 

3 . 3 Listings 

Subroutine SOLVE is self-explanatory and needs no elaboration, 
except as to the manner in which the equations must be stored. The 
blocks are stored sequentially by rows, as shown below. At the end 
of each row, a block of load vectors is present: 
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K 11 K 12 K 13 K 14 K lm R 1 

K 22 K 23 K 24 K 25 K 2 m R 2 



K II K I,I+1 



Rj (37) 



K . - K R 

n-l,n-l n-l,n n-1 



The symbols m and n are the number of blocks per row and column, 
respectively. 

This method of storage leaves a few blocks of unused storage. 
However, the last row of blocks is needed for temporary storage during 
the elimination. Three service subroutines are needed in SOLVE. The 
multiplication and inversion operations are standard and need no words 
of explanation. The subroutine WRDISK/RDDISK, however, is hardware 
dependent and uses FORTRAN commands particular to IBM/360 installa- 
tions (Reference 14). This feature, although undesirable, is not a very 
serious drawback because every installation using fast random access 
devices has specialized commands to use these facilities. For other 
installations, this subroutine would have to be rewritten to conform 
with the available hardware. The listings are self-explanatory. 
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3 . 4 Timing 



The use of the program listed in Appendix 4 poses a serious practi- 
cal problem. The solution of large systems of equations requires a 
large amount of computer time. The most important factors affecting 
the required time in the system proposed are: the number of calls to 

the WRDISK/RDDISK routine, and the size of the blocks (N ) . In the 

s 

program, the number of calls to WRDISK/RDDISK has been reduced as 
much as we could. In actual tests, the object codes were optimized, 
with respect to execution time, by the compilers that produced these 
object codes. We used the H compiler of release 18 of the IBM/360 
operating system. A system of 500 equations with a half band vyidth 
of 151 was solved with three different block sizes, (25, 50, 100). 

The execution times were (510, 402, 297) seconds, respectively. On 
the basis of this test, it seems that the block-size should be as large 
as the computer will allow, in order to reduce the execution time. For 
a block-size of 50, the following empirical formula gives an estimate 
of the required time: 

T = 1.857 nm 2 + 0.102 n 2 m - 0.322 n 2 + 4.682 nm- 5.230 m 2 (38) 
The formula was obtained with the tests shown in Table 1. For the 
first test n = 4, m = 3, it was possible to obtain an in-core solution 
with the subroutine shown in Appendix 3. The execution time was 61.3 
seconds for the in-core solution, as compared to 7 5.7 seconds for the 
block solver. This is indeed a very small penalty if we also compare the 
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storage requirements of the two solvers; 19 6K bytes for the in-core 
solver and 122K bytes for the block solver. 

In a recent study. Reference (15), various equation solvers were 
compared. Using a sliding block technique similar to that mentioned 
in Reference (8), a system of 1080 equations with a half-band width 
of 357 was solved in 1365.1 seconds. The machine used was a CDC 
6600 and the method used required 73728 words of core storage. 

Test Number 5 in Table 1 shows that 1000 equations with a half -band 
width of 351 required 2989.8 seconds and used 15250 double words of 
core storage. Again, this is a comparison that seems very favorable 
to the block solver, if we keep in mind that the CDC6600 is a much 
faster machine than the IBM/360/67. Average execution times for 
division, multiplication and addition are: (2.9, 1.0, 0.4) micro- 
seconds for the CDC machine and (13.9, 10.6, 8.0) for the IBM machine. 
3 . 5 Accuracy 

The most critical operation, with respect to round-off, is the 

inversion of the diagonal blocks. If the block is singular, the program 

will print a diagnostic to this effect. However, this is not likely to 

happen very often. A much more serious difficulty comes from blocks 

that are nearly singular and this must also be detected. For this 

reason, the inversion routine will calculate an estimate of the round-off 

-12 

value of zero as the trace of the block multiplied by 10 and report 
that the block is nearly singular, if any diagonal term becomes smaller 
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TABLE 1 TIMING TESTS 



n 

(Number of 
blocks per 
column) 


m 

(Number of 
blocks per 
row) 


N 

(Number of 
equations) 


M 

(Half-band 

width) 


T 

(Time in seconds) 


4 


3 


200 


101 


75.7 


100 


3 


5000 


101 


2977.4 


10 


5 


500 


201 


586.5 


50 


5 


2500 


201 


3834.8 


20 


8 


1000 


351 


2989.8 



* 

The blocksize used was 50 and the total storage requirement of the 
test program was 122K bytes or 15250 double words. 



than this value during inversion. Besides these diagnostics, the 
program will also print a condition number for each diagonal block. 

This condition number is obtained as the product of the maximum column 
sum of the block before and after inversion, following Wilkinson's 
perturbation theory. With this condition number, an estimate of the 
relative error in the answer can be obtained from: 



E(%) £ c 




100 



( 39 ) 



where c is the maximum value of all condition numbers of the system, 

||p|| is the sum of the absolute values of the estimated errors in the 
loads and ||r|| is the sum of the absolute values of the loads. 

3 . 6 Modifications 

The program is presented in its simplest form; no facilities for • 
iteration and multiple-load vectors are provided for. Users with particu- 
lar needs should have no difficulty modifying the program listed, in 
order to take advantage of hardware available in their computing 
facilities . 

3 . 7 Conclusion 

The program presented here provides a very simple answer to a 
most difficult practical difficulty. When a problem is too big for a 
band or a frontal solver, the block solver will still give a solution, 
even with a very small computer, provided one is ready to pay a 
penalty in execution time . 
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Chapter IV. Conclusions 



The publication of this report was delayed in order to verify the 
integrity of the service routines that it contains. This has now been 
completed and it can be reported that the routines perform as expected. 

4 . 1 Stress-strain analysis codes 

Several computer codes for mesh generation and consistent load 
vectors computation, as well as two basic stress strain programs for 
two and three dimensional problems, have been completed and tested 
(16), (17) and (18). From that work it may be concluded that these ele- 
ments will shed significant new light on the stress analysis of rocket 
motors . 

4.2 Reduced Gaussian Integration 

In both reference (17) and (18), it is reported that the Gaussian 
integration mentioned in article 1.9 leads to improved results when 
only 2 and 3 Gauss ordinates are uses, respectively, in the generation 
of the 20 and 32 nodal point elements. This reduced integration also 
leads to significant economies in the total solution times. 

4 . 3 Surface stresses 

In some problems it was observed that the stress and strain values 
obtained at the nodes were in error on the external surfaces of the ob- 
jects analysed. Internal nodes, however, gave the best results. This 
matter deserves further investigations . 

4 . 4 Future work 

The program^developed during the course of this investigation will 
be refined and released as soon as they are ready. In the meantime, 
potential users can refer to reference (17) for a first release of the three- 
dimensional programs . 
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Appendix I. 



Code for the generation of two-dimensional elements 
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SUBROUTINE C LR *♦ ( S T K ♦ A K , 1} t N ) C4P0CC10 

CUF4 IS A NUMERICAL I NTFGR AT I PN SLBRCUTTNF THAT FORMS THF STIFFNESS 
MATRIX FOR A QUADRATIC THREE DIMENSIONAL ISOPARAMETRIC ELEMENT. 

FOUR GAUSS ORDINATES ARF USED. 
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Appendix 3 . Listing for an in-core band solver 
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SUBROUTINE SYMSOLI A , B ,C ♦ NN , MM ) 
IMPLICIT REAL*3( A-H.O-Z ) 
DIMENSION A(NN,MM) , B(NN) ,C(MM) 
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Appendix 4. Listing for the Block Solver 
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